In [24]:
import xarray as xr
import polars as pl
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import Point
import numpy as np
import pandas as pd
import calendar
from shapely.geometry import LineString
from matplotlib.cm import ScalarMappable
from matplotlib.colors import Normalize, LinearSegmentedColormap
# taken from https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r01/access/netcdf/IBTrACS.ALL.v04r01.nc on 11/4/2025
filepath = "../.data/IBTrACS.ALL.v04r01.nc"

ds = xr.open_dataset(filepath)

country1 = 'Jamaica'
country1_alt = ""
distance = 50
MILES_TO_METERS = 1609.34
BUFFER_DISTANCE_MILES = distance
BUFFER_DISTANCE_METERS = BUFFER_DISTANCE_MILES * MILES_TO_METERS
MPH_TO_KNOTS = 1.15078
country_string = country1
if country1_alt == "":
    country_string_alt = country1
else:
    country_string_alt = country1_alt
country_URL = "../.data/ne_50m_admin_0_countries.zip" # from https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/ne_50m_admin_0_countries.zip
country_gdf = gpd.read_file(country_URL)
country_shape = country_gdf[country_gdf['ADMIN'] == country_string_alt]
country_shape_proj = country_shape.to_crs("EPSG:3857")
country_buffer_proj = country_shape_proj.geometry.buffer(BUFFER_DISTANCE_METERS)
country_buffer_gdf = gpd.GeoDataFrame(geometry=country_buffer_proj, crs="EPSG:3857")
country_buffer_gdf = country_buffer_gdf.to_crs("EPSG:4326")
country_geometry = gpd.GeoDataFrame(
    geometry=[country_buffer_gdf.union_all()], 
    crs="EPSG:4326"
)
country_shape_original = country_gdf[country_gdf['ADMIN'] == country_string_alt].to_crs("EPSG:4326")
country_geometry_original = gpd.GeoDataFrame(
    geometry=[country_shape_original.union_all()], 
    crs=country_shape_original.crs
)
# country_geometry = gpd.GeoDataFrame(geometry=[country_shape.union_all()], crs=country_shape.crs)
track_data = ds[['sid', 'time', 'lat', 'lon', 'usa_wind', 'usa_sshs', 'storm_speed']]
track_data['usa_wind'] = track_data['usa_wind'] * MPH_TO_KNOTS
track_data['storm_speed'] = track_data['storm_speed'] * MPH_TO_KNOTS
track_df = track_data.to_dataframe().reset_index()
track_df = track_df.dropna(subset=['lon', 'lat'])
geometry = [
    Point(lon, lat) 
    for lon, lat in zip(track_df['lon'], track_df['lat'])
]
TARGET_CRS = "EPSG:4326"
storms_gdf = gpd.GeoDataFrame(track_df, geometry=geometry, crs=TARGET_CRS)
storms_over_country_gdf = gpd.sjoin(
    storms_gdf, 
    country_geometry, 
    how="inner", 
    predicate="within"
)
storms_gdf['is_inside_country'] = storms_gdf.index.isin(storms_over_country_gdf.index)
storms_pl_full = pl.from_pandas(storms_gdf.drop(columns=['geometry']))
storms_pl_full = storms_pl_full.sort(["storm", "date_time"])
storms_inside_area = storms_pl_full.filter(
    pl.col('is_inside_country')
)
max_intensity_within_area = (
    storms_inside_area
    .filter(pl.col('usa_wind').is_not_null())
    .sort(['storm', 'usa_wind'], descending=[False, True]) 
    .unique(subset=['storm'], keep='first') 
    .select('storm', 'date_time', 'lon', 'lat', 'usa_wind', 'usa_sshs', 'time', 'sid', 'storm_speed')
)
max_intensity_within_area_pandas = max_intensity_within_area.to_pandas()
geometry = [
    Point(lon, lat) 
    for lon, lat in zip(max_intensity_within_area_pandas['lon'], max_intensity_within_area_pandas['lat'])
]
country_cyclones_pandas = max_intensity_within_area_pandas
geometry = [
    Point(lon, lat)
    for lon, lat in zip(country_cyclones_pandas['lon'],
                        country_cyclones_pandas['lat'])
]
country_cyclones_gdf = gpd.GeoDataFrame(
    country_cyclones_pandas,
    geometry=geometry,
    crs="EPSG:4326"
)

fig, ax = plt.subplots(1, 1, figsize=(10, 8))

colors = ["#192D66", "#FF00FF"]  # Dark blue to Magenta
custom_cmap = LinearSegmentedColormap.from_list("dark_blue_magenta", colors)

country_geometry_original.plot(
    ax=ax, 
    color='lightgreen',              # Fill color for the land
    edgecolor='darkgreen',           # Border color
    linewidth=1.5,
    zorder=1,
    alpha=0.3
)
country_cyclones_gdf.plot(
    ax=ax,
    marker='o',
    column=country_cyclones_gdf['usa_wind'],
    cmap=custom_cmap,                   
    markersize=50,                  
    legend=True,
    legend_kwds={'label': "Wind Speed (mph)", 'orientation': "horizontal"},
    zorder=3
)
minx, miny, maxx, maxy = country_cyclones_gdf.total_bounds
if country_string == "United States":
    ax.set_xlim(-115, -60)
    ax.set_ylim(16, 45)
else:
    ax.set_xlim(minx - 0.5, maxx + 0.5)
    ax.set_ylim(miny - 0.5, maxy + 0.5)

ax.set_title(f"Storm Track Intensity Over {country_string} (Color by Max Wind Speed)", fontsize=14)
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")

plt.tight_layout()
plt.show()
No description has been provided for this image
In [25]:
mean_wind = country_cyclones_pandas['usa_wind'].mean()
std_wind = country_cyclones_pandas['usa_wind'].std()
ci_upper = mean_wind + 1.96 * std_wind
ci_lower = mean_wind - 1.96 * std_wind

fig, ax = plt.subplots(1, 1, figsize=(10, 8))

plt.scatter(country_cyclones_pandas['time'], country_cyclones_pandas['usa_wind'])

ax.set_title(f"Tropical Cyclone Strength Over Time for {country_string} (mean and 95% CI shown)", fontsize=14)
ax.set_xlabel("Time")
ax.set_ylabel("Wind Speed (mph)")
ax.axhline(
    mean_wind, 
    color="#FF00FF", 
    linestyle='-', 
    linewidth=2.5, 
    label=f'Overall Mean Wind Speed ({mean_wind:.1f} mph)'
)
ax.axhline(
    ci_upper, 
    color="#FF00FF", 
    linestyle='--', 
    linewidth=1.5, 
    label=f'95% CI Upper Bound ({ci_upper:.1f} mph)'
)
ax.axhline(
    ci_lower, 
    color="#FF00FF", 
    linestyle='--', 
    linewidth=1.5,
    label=f'95% CI Lower Bound ({ci_lower:.1f} mph)'
)

plt.tight_layout()
plt.show()
No description has been provided for this image
In [26]:
jamaica_wind_speeds = country_cyclones_pandas[['sid','usa_wind']].rename(columns={'usa_wind': 'max_storm_wind'})
jamaica_wind_speeds['sid'] = jamaica_wind_speeds['sid'].str.decode('utf-8')
jamaica_wind_speeds['sid'] = pd.Series(jamaica_wind_speeds['sid'].tolist())
relevant_storms = country_cyclones_pandas[country_cyclones_pandas['usa_wind'] > 75]['sid'].unique()
# target_sids_list = [sid.decode('utf-8') for sid in relevant_storms]
mask = ds['sid'].isin(relevant_storms)
ds_filtered = ds.where(mask, drop=True)
track_data_vars = ['sid', 'time', 'lat', 'lon', 'usa_wind']
track_df = ds_filtered[track_data_vars].to_dataframe().reset_index()
track_df['sid'] = track_df['sid'].str.decode('utf-8')
track_df['sid'] = pd.Series(track_df['sid'].tolist())
track_df = track_df.sort_values(by=['sid', 'time']).reset_index(drop=True)
# track_df = track_df.rename(columns={'sid': 'storm', 'time': 'date_time'})
track_df = track_df.dropna(subset=['lon', 'lat', 'usa_wind'])
geometry = [
    Point(lon, lat) 
    for lon, lat in zip(track_df['lon'], track_df['lat'])
]
TARGET_CRS = "EPSG:4326"
storms_gdf = gpd.GeoDataFrame(
    track_df, 
    geometry=geometry, 
    crs=TARGET_CRS
)
# map_shape = country_gdf
# map_proj = map_shape.to_crs("EPSG:3857")
# map_gdf = gpd.GeoDataFrame(geometry=map_proj, crs="EPSG:3857")
# map_gdf = map_gdf.to_crs("EPSG:4326")
# map_geometry = gpd.GeoDataFrame(
#     geometry=[map_gdf.union_all()], 
#     crs="EPSG:4326"
# )
map_shape = country_gdf.to_crs("EPSG:4326")
map_geometry = gpd.GeoDataFrame(
    geometry=[map_shape.union_all()], 
    crs=map_shape.crs
)
storms_gdf = storms_gdf.sort_values(by=['sid', 'time']) 

# max_wind_per_storm = storms_gdf.groupby('sid')['usa_wind'].max().rename('max_storm_wind').reset_index()
lines_series = (
    storms_gdf.groupby('sid')['geometry']
    .agg(lambda x: LineString(x.tolist()) if len(x.tolist()) > 1 else None)
    .rename('geometry')
)

storms_lines_gdf = gpd.GeoDataFrame(lines_series).reset_index()
storms_lines_gdf = storms_lines_gdf.set_geometry('geometry', crs=storms_gdf.crs)
storms_lines_gdf = storms_lines_gdf.merge(jamaica_wind_speeds, on='sid', how='left')

fig, ax = plt.subplots(1, 1, figsize=(12, 10))

map_geometry.plot(ax=ax, color='lightgreen', edgecolor='darkgreen', linewidth=1.5, zorder=1, alpha=0.3)
ax.set_xlim(-100, -50)
ax.set_ylim(10, 40)
ax.set_aspect('equal')

colors = ["#192D66", "#FF00FF"]  # Dark blue to Magenta
custom_cmap = LinearSegmentedColormap.from_list("dark_blue_magenta", colors)

storms_lines_gdf.plot(
    ax=ax,
    column='max_storm_wind', # Color the line by its max wind speed (the merged column)
    cmap=custom_cmap,
    legend=False,
    legend_kwds={'label': "Max Storm Wind Speed (mph)", 'orientation': "horizontal"},
    linewidth=2.0,
    linestyle='-',
    zorder=2
)

country_geometry_original.plot(
    ax=ax, 
    color='lightgreen',              # Fill color for the land
    edgecolor='darkgreen',           # Border color
    linewidth=1.5,
    zorder=5,
    alpha=0.9
)

strongest = storms_lines_gdf.sort_values(by='max_storm_wind', ascending=False).head(1)

strongest.plot(
    ax=ax,
    color="#FF00FF",
    legend=False,
    legend_kwds={'label': "Max Storm Wind Speed (mph)", 'orientation': "horizontal"},
    linewidth=5.0,
    linestyle='-',
    zorder=3
)
min_val = storms_lines_gdf['max_storm_wind'].min()
max_val = storms_lines_gdf['max_storm_wind'].max()
sm = ScalarMappable(cmap=custom_cmap, norm=Normalize(vmin=min_val, vmax=max_val))
sm.set_array([])
cbar = fig.colorbar(sm, ax=ax, orientation='horizontal', fraction=0.05, pad=0.1)
cbar.set_label('Max Storm Wind Speed (mph)', fontsize=12)

ax.set_title("Historical Hurricane Tracks with over 75 mph Wind (Line Color by Max Wind)", fontsize=16)
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
plt.show()
No description has been provided for this image
In [27]:
month_analysis = pl.DataFrame(country_cyclones_pandas)
month_analysis = (
    month_analysis
    .with_columns(
        pl.col('time').dt.month().alias('month')
    )
)

months = month_analysis.group_by('month').agg(pl.len().alias('storm_count')).sort('month').to_pandas()


fig, ax = plt.subplots(1, 1, figsize=(10, 6))

bars = ax.bar(
    x=months['month'], 
    height=months['storm_count'],
    color='#1f77b4',       # Standard Matplotlib blue (can be customized)
    edgecolor='black',
    linewidth=0.5
)

month_names = [calendar.month_abbr[i] for i in months['month']]
ax.set_xticks(months['month'])
ax.set_xticklabels(month_names)

for bar in bars:
    yval = bar.get_height()
    ax.text(
        bar.get_x() + bar.get_width()/2, 
        yval + 0.5,
        int(yval), 
        ha='center', 
        va='bottom'
    )

ax.set_title(
    'Jamaica Tropical Cyclone Frequency by Month (Storms with over 75 mph Wind)', 
    fontsize=16, 
    fontweight='bold'
)
ax.set_xlabel('Month', fontsize=12)
ax.set_ylabel('Number of Storms', fontsize=12)

ax.set_ylim(0, months['storm_count'].max() * 1.1)

ax.grid(axis='y', linestyle='--', alpha=0.6)

ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

plt.tight_layout()
plt.show()
No description has been provided for this image
In [28]:
month_analysis = pl.DataFrame(country_cyclones_pandas)
month_analysis = (
    month_analysis
    .filter(pl.col('usa_wind') > 100)
    .with_columns(
        pl.col('time').dt.month().alias('month')
    )
)

months = month_analysis.group_by('month').agg(pl.len().alias('storm_count')).sort('month').to_pandas()


fig, ax = plt.subplots(1, 1, figsize=(10, 6))

bars = ax.bar(
    x=months['month'], 
    height=months['storm_count'],
    color='#1f77b4',       # Standard Matplotlib blue (can be customized)
    edgecolor='black',
    linewidth=0.5
)

month_names = [calendar.month_abbr[i] for i in months['month']]
ax.set_xticks(months['month'])
ax.set_xticklabels(month_names)

for bar in bars:
    yval = bar.get_height()
    ax.text(
        bar.get_x() + bar.get_width()/2, 
        yval + 0.2,
        int(yval), 
        ha='center', 
        va='bottom'
    )

ax.set_title(
    'Jamaica Tropical Cyclone Frequency by Month (Storms with over 100 mph Wind)', 
    fontsize=16, 
    fontweight='bold'
)
ax.set_xlabel('Month', fontsize=12)
ax.set_ylabel('Number of Storms', fontsize=12)

ax.set_ylim(0, months['storm_count'].max() * 1.1)

ax.grid(axis='y', linestyle='--', alpha=0.6)

ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

plt.tight_layout()
plt.show()
No description has been provided for this image
In [29]:
enso_index = pl.read_csv("../.data/ENSO_index.csv") #from https://www.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php
enso_long = (
    enso_index.unpivot(
        index=['Year'],
        on=enso_index.columns[1:],
        variable_name='past_3_months',
        value_name='ENSO_Index'
    )
    .rename({'Year':'year'})
)
ENSO_map = pl.DataFrame({
    "past_3_months": [
        'OND',
        'NDJ',
        'DJF',
        'JFM',
        'FMA',
        'MAM',
        'AMJ',
        'MJJ',
        'JJA',
        'JAS',
        'ASO',
        'SON',],
    "month": [1,2,3,4,5,6,7,8,9,10,11,12]
})
enso_analysis = (
    pl.DataFrame(country_cyclones_pandas)
    .with_columns(
        pl.col('time').dt.month().alias('month'),
        pl.col('time').dt.year().alias('year')
    )
    .with_columns(
        pl.when(pl.col('month').is_in([1,2])) # (unlikely to occur)
        .then(pl.col('year')-1) # this is because DJF is in the current year, but OND and NDJ are not.
        .otherwise('year')
        .alias('year')
    )
    .join(ENSO_map, on='month', how='left')
    .join(enso_long, on=['year','past_3_months'], how='left')
    .with_columns(
        pl.when(pl.col('month').is_in([1,2]))
        .then(pl.col('year') + 1) # changing back to actual year after merge
        .otherwise('year')
        .alias('year')
    )
)

fig, ax = plt.subplots(1, 1, figsize=(10, 6))

ax.scatter(
    enso_analysis['ENSO_Index'],
    enso_analysis['usa_wind'],
    color='#1f77b4',
    alpha=0.7,
    edgecolors='black',
    s=50,
)


ax.set_title(
    'Relationship Between Prior 3 month ENSO Index and Tropical Cyclone Wind Speed (mph)',
    fontsize=14,
    fontweight='bold'
)
ax.set_xlabel('ENSO Index', fontsize=12)
ax.set_ylabel('Wind Speed (mph)', fontsize=12)

ax.grid(axis='both', linestyle='--', alpha=0.6)

plt.tight_layout()
plt.show()
No description has been provided for this image
In [30]:
# API Request for ERA5 data:
# import cdsapi

# dataset = "reanalysis-era5-single-levels"
# request = {
#     "product_type": ["reanalysis"],
#     "variable": ["sea_surface_temperature"],
#     "year": [
#         "1940", "1941", "1942",
#         "1943", "1944", "1945",
#         "1946", "1947", "1948",
#         "1949", "1950", "1951",
#         "1952", "1953", "1954",
#         "1955", "1956", "1957",
#         "1958", "1959", "1960",
#         "1961", "1962", "1963",
#         "1964", "1965", "1966",
#         "1967", "1968", "1969",
#         "1970", "1971", "1972",
#         "1973", "1974", "1975",
#         "1976", "1977", "1978",
#         "1979", "1980", "1981",
#         "1982", "1983", "1984",
#         "1985", "1986", "1987",
#         "1988", "1989", "1990",
#         "1991", "1992", "1993",
#         "1994", "1995", "1996",
#         "1997", "1998", "1999",
#         "2000", "2001", "2002",
#         "2003", "2004", "2005",
#         "2006", "2007", "2008",
#         "2009", "2010", "2011",
#         "2012", "2013", "2014",
#         "2015", "2016", "2017",
#         "2018", "2019", "2020",
#         "2021", "2022", "2023",
#         "2024", "2025"
#     ],
#     "month": [
#         "01", "02", "03",
#         "04", "05", "06",
#         "07", "08", "09",
#         "10", "11", "12"
#     ],
#     "day": ["15"],
#     "time": ["00:00"],
#     "data_format": "netcdf",
#     "download_format": "zip",
#     "area": [17.5, -77, 16.5, -76]
# }

# client = cdsapi.Client()
# client.retrieve(dataset, request).download()

sst_ds = xr.open_dataset("../.data/sst_caribbean.nc")

sst_pl = (
    pl.DataFrame(
        sst_ds
        .mean(dim=['latitude', 'longitude'])
        .to_pandas()
        .drop(columns=['expver','number'])
        .reset_index()
    )
    .with_columns(
        pl.col('valid_time').dt.month().alias('month'),
        pl.col('valid_time').dt.year().alias('year'),
    )
    .drop('valid_time')
)

sst_analysis = (
    pl.DataFrame(country_cyclones_pandas)
    .with_columns(
        pl.col('time').dt.month().alias('month'),
        pl.col('time').dt.year().alias('year')
    )
    .join(sst_pl, on=['month', 'year'], how='left')
    .with_columns(
        (pl.col('sst') - 273.15).alias('sst')
    )   
)

fig, ax = plt.subplots(1, 1, figsize=(10, 6))

ax.scatter(
    sst_analysis['sst'],
    sst_analysis['usa_wind'],
    color='#1f77b4',
    alpha=0.7,
    edgecolors='black',
    s=50,
)

ax.set_title(
    'Relationship Between Sea Surface Temperature (°C) and Tropical Cyclone Wind Speed (mph)',
    fontsize=14,
    fontweight='bold'
)
ax.set_xlabel('SST (°C)', fontsize=12)
ax.set_ylabel('Wind Speed (mph)', fontsize=12)

ax.grid(axis='both', linestyle='--', alpha=0.6)

plt.tight_layout()
plt.show()
No description has been provided for this image
In [31]:
fig, ax = plt.subplots(1, 1, figsize=(10, 6))

ax.scatter(
    sst_analysis['storm_speed'],
    sst_analysis['usa_wind'],
    color='#1f77b4',
    alpha=0.7,
    edgecolors='black',
    s=50,
)

ax.set_title(
    'Relationship Between Storm Speed (mph) and Tropical Cyclone Wind Speed (mph)',
    fontsize=14,
    fontweight='bold'
)
ax.set_xlabel('Storm Speed (mph)', fontsize=12)
ax.set_ylabel('Wind Speed (mph)', fontsize=12)

ax.grid(axis='both', linestyle='--', alpha=0.6)

plt.tight_layout()
plt.show()
No description has been provided for this image
In [32]:
# Ivan
pl.DataFrame(country_cyclones_pandas).filter((pl.col('time').dt.year() == 2004) & (pl.col('time').dt.month() == 9))
Out[32]:
shape: (1, 9)
stormdate_timelonlatusa_windusa_sshstimesidstorm_speed
i64i64f32f32f32f32datetime[ns]binaryf32
1134167-76.517.299999155.3553014.02004-09-11 00:00:00.000040448b"2004247N10332"11.507799
In [ ]:
# Dean
pl.DataFrame(country_cyclones_pandas).filter((pl.col('time').dt.year() == 2007) & (pl.col('time').dt.month() == 8))
Out[ ]:
shape: (1, 9)
stormdate_timelonlatusa_windusa_sshstimesidstorm_speed
i64i64f32f32f32f32datetime[ns]binaryf32
1163953-76.90000217.299999143.8474884.02007-08-19 21:00:00.000039936b"2007225N12331"20.714039
In [34]:
# Beryl
pl.DataFrame(country_cyclones_pandas).filter((pl.col('time').dt.year() == 2024) & (pl.col('time').dt.month() == 7))
Out[34]:
shape: (1, 9)
stormdate_timelonlatusa_windusa_sshstimesidstorm_speed
i64i64f32f32f32f32datetime[ns]binaryf32
1335244-76.80000317.299999138.0935974.02024-07-03 18:00:00.000039936b"2024181N09320"17.2617
In [35]:
# Gilbert
pl.DataFrame(country_cyclones_pandas).filter((pl.col('time').dt.year() == 1988) & (pl.col('time').dt.month() == 9))
Out[35]:
shape: (1, 9)
stormdate_timelonlatusa_windusa_sshstimesidstorm_speed
i64i64f32f32f32f32datetime[ns]binaryf32
959632-76.80000317.799999132.3396914.01988-09-12 17:00:00.000053760b"1988253N12306"13.80936
In [36]:
# Charlie
pl.DataFrame(country_cyclones_pandas).filter((pl.col('time').dt.year() == 1951) & (pl.col('time').dt.month() == 8))
Out[36]:
shape: (1, 9)
stormdate_timelonlatusa_windusa_sshstimesidstorm_speed
i64i64f32f32f32f32datetime[ns]binaryf32
531248-76.19999717.6126.5857933.01951-08-18 00:00:00.000040448b"1951224N12316"17.2617
In [37]:
pl.DataFrame(country_cyclones_pandas).sort('usa_wind')
Out[37]:
shape: (66, 9)
stormdate_timelonlatusa_windusa_sshstimesidstorm_speed
i64i64f32f32f32f32datetime[ns]binaryf32
880156-77.017.417.2617-3.01981-08-14 06:00:00.000040448b"1981219N11334"14.960139
1132124-76.19999717.20000128.769499-3.02004-08-06 12:00:00.000040448b"2004217N13306"13.80936
114750-78.517.628.769499-1.02005-10-15 18:00:00.000039936b"2005289N18282"3.45234
26936-76.80000317.20000134.523399-1.01916-10-10 00:00:00.000040192b"1916283N19284"6.90468
95123-78.90000217.934.523399-1.01987-11-01 03:00:00.000040448b"1987305N17283"13.80936
………………………
959632-76.80000317.799999132.3396914.01988-09-12 17:00:00.000053760b"1988253N12306"13.80936
1335244-76.80000317.299999138.0935974.02024-07-03 18:00:00.000039936b"2024181N09320"17.2617
1163953-76.90000217.299999143.8474884.02007-08-19 21:00:00.000039936b"2007225N12331"20.714039
1134167-76.517.299999155.3553014.02004-09-11 00:00:00.000040448b"2004247N10332"11.507799
1352289-78.017.9184.1247865.02025-10-28 15:00:00.000039936b"2025291N11319"8.05546
In [ ]: